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A COMPUTER PROGRAM FOR RAY TRACE 


THROUGH SPACECRAFT WINDOWS 
Kenneth C. White 
Ames Research Center 


SUMMARY 


A FORTRAN IV computer program has been written specifically for tracing 
rays through spacecraft windows. The program computes the window-induced 
angular deviation of the rays from their original path. Any type of window 
may be considered with no restriction on size, shape, material, or number of 
panes. The necessary equations are written in vector matrix form for mathe- 
matical convenience and ease in computation. The program requires mathemat- 
ical models of each window surface. The numerical procedures used are 
described and a test case is presented. 


INTRODUCTION 


On-board navigation systems that derive their basic inputs from optical 
measurements made through a spacecraft window have been considered for manned 
space flight. The optical measurements of interest are the angle between two 
celestial bodies, such as a star and a planet, and the angle between a star 
and a spacecraft. These angular navigation measurements are subject to 
window-induced errors that result from the angular deviation or bending of the 
lines of sight as they pass through the window. In this report the terms 
"ray" and "line of sight" are used interchangeably. It should be noted that 
ray tracing is done in the opposite direction of the normal ray tracing. 

The purpose of this program is to compute the angular line-of-sight 
deviations induced by the spacecraft windows. This is accomplished by tracing 
the lines of sight through the window by geometric ray tracing techniques. 

Many ray trace programs are available but they are primarily oriented toward 
lens design and usually require rotationally symmetric optical systems. A 
FORTRAN IV computer program has therefore been written which emphasizes the 
computation of the angular deviations of the rays and is not limited to 
rotationally symmetric optical systems. 

The input to the program consists of mathematical models of each window 
surface, the orientation parameters of the lines of sight to be traced through 
the window, and the various parameters of the window environment. The output 
consists of the orientation parameters of the line of sight after it has 
passed through the window and the angular differences between the original 
line of sight and the deviated line of sight. 



The basic equations are presented and the program is described fully in 
this report. This description includes the program listing, program usage, 
flow charts, and a sample case. 

An application of the program and a discussion of analytical methods 
which utilize the program are given in reference 1. 


PROGRAM DESCRIPTION 


The presence of a spacecraft window in the path of rays of light causes 
these rays to bend or to deviate from their original path. This program is 
designed to determine the magnitude and direction of the angular deviations of 
the light rays. In order to determine the deviations, the ray must be traced 
through the window and a comparison made between the entering and exiting rays. 
The fundamental law of ray tracing is Snell's law of refraction, 

ni sin 0^ = n 2 sin 0 2 


where 

nj the index of refraction of medium 1 

01 the angle the incident ray makes with the normal to the interface 

between mediums 1 and 2 

n 2 the index of refraction of medium 2 

0 2 the angle the refracted ray in medium 2 makes with the normal 

Snell's law is used to trace a ray through several mediums by successive 
applications of the law at each interface between mediums such as air and 
glass. This is accomplished by letting the refracted ray become the incident 
ray for the next application of the law. In order to apply Snell's law, the 
normals to the interfaces between mediums at the points where the ray inter- 
sects the interface must be known. This means the shape of the interface must 
be known. An iteration scheme is also required to determine the points of 
intersection of the ray with the interfaces. 


BASIC EQUATIONS 


For mathematical convenience and ease in computation, Snell's law is 
expressed in vector form as follows: 
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where I and R are the unit vectors in the direction of the incident and 
refracted rays. The unit vector N is in the direction of the normal to the 
interface between mediums 1 and 2 at the point of intersection of the incident 
ray; ni and n 2 are as previously defined. Equation (1) is derived in 
reference 2. 


Coordinate System 

The coordinate system is illustrated in figure 1. The innermost surface 
of the spacecraft window is assumed to lie in the xy-plane. The positive 

+y z-axis is toward the out- 

f side of the spacecraft. 

The incident and refracted 
rays and their orientation 
'3*4 \||. / angles are also illustrated 

X. in the figure. The vector 

8 . i \ I is the unit vector in 

yg * J^f P the direction of the inci- 

&/\/ /v : 3r r|. dent ray. The azimuth 

angle, a^, is defined as 

\ \C| X ::y| the angle in the xy-plane 

\ \\ M/ \ m between the positive 

\\ \ x-axis and the projection, 

V V*r \<*i Ip, of I onto the xy- 

\ \ \ plane. The elevation angle, 

\ \ 1188 1 . , is the angle between 



and I T 


The incidence 


emerges from the outer-most 

Figure 1.- Ray trace coordinate system. window Surface. The azi- 

muth angle, a r , is the 

angle in the xy-plane between the positive x-axis and the proj ection, ^Rp , of 
R onto the xy-plane. The elevation angle, 6 r , is the angle between R and 
Rp. The azimuth angles and a r are measured from the positive x-axis 
toward the positive y-axis and vary from 0° to 360° x The elevation angles. 


Figure 1.- Ray trace coordinate system. 


6 q an d 6 ) 
positive 


are measured from the projected vectors 
z-axis and vary from 0° to 90°. 


and R-n toward the 


The unit vectors I and R are given by the following equations: 

^ ✓>. A A 

I = cos cos i + cos 6^ sin ctp j + sin 6^ k 


R = cos 6 r cos a r i + cos 6 sin a r j + sin 6 r k 


where i, j, and k are unit vectors in the directions of the x, y, and z 
axes, respectively. 
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Two angular ray deviations Aa and A6 are also defined in figure 1. 

The out-of-plane deviation is defined as Aa = (a^ - a r ) . The in-plane devia- 
tion is defined as A6 = (6p - 6 r ) . These angular deviations are the 
important results of the ray trace. 


Window Surface Mathematical Models 

The window surface shapes are described in terms of their deviations from 
reference planes set parallel to the xy-plane. The first reference plane is 
coincident with the xy-plane. There is a reference plane for each window 
surface. The window surface deviation from the reference plane AZ is given 
by a fourth degree mixed polynomial in x and y as follows: 


AZ = m n x 4 y 4 + m 12 x 3 y 4 + m 13 x 2 y 4 + • • • + m 45 Y + m 55 


(4) 


or in shorthand form: 


_T _ 
AZ = Y MX 


where 


Y t = [y 4 y 3 y 2 y 1] 



x 

1 


mu 

m l2 

mi 3 

m 14 

m 15 

m 21 

m 22 

m 2 3 

m 24 

m 2 5 

m 3 1 

m 32 

m 33 

m 34 

m 35 | 

m 41 

m 42 

m 43 

m 44 

m 45 

L?*51 

m 52 

m 53 

m 54 

m 5 5_ 


(5) 


The elements of the matrix M are program inputs. The normal to the surface 
at a given point is obtained by evaluating the gradient of the surface at that 
point. If equation (5) is written in the form 

F(x, y, z) = AZ - Y T MX (6) 
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then the unit vector normal to the surface F is given by 



where 


( 8 ) 


Equations (4) through (8) are derived in reference 2. 

Equations (1) through (8) together with the iteration scheme described in 
appendix A comprise the necessary elements to compute a ray trace. 

Description of a Ray Trace 

A block diagram illustrating the logical sequencing of the ray trace 
program is given in figure 2. The orientation angles of the incident ray, the 
coordinates of the reference point, the mathematical models of the window sur- 
faces, and the window and environmental parameters are input. The components 
of the unit vector I in the direction of the incident^ ray are computed by 
means of equation (2). The components of the vector E from the origin to 
the reference point are computed. Next the point of intersection of I with 
the first reference plane is computed using equations in appendix A. Then the 
point of intersection of I with the first window surface is calculated using 
the iteration scheme described in appendix A. The iteration is done in 
subroutine ITER. After the point of intersection is found, the unit vector N 
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Figure 2.- Computer program schematic. 


normal to the window surface at the point is computed by means of equations (6) , 
(7), and (8). This is accomplished in subroutine NORMAL. The refracted ray 
in the next medium is computed by equation (1). This is done in subroutine 
REFRC. The number of surfaces traversed by the ray are tested. If the ray 
has passed through the last surface, the refracted ray is compared with the 
original incident ray to determine the angular deviations. If the ray has not 
passed through all the surfaces, the refracted ray becomes the incident ray 
and the ray trace is continued until the ray has completely passed through the 
window. 


ROUTINE DESCRIPTIONS 


The main program is designated MS2500 and is the main processing unit for 
the ray trace. Equation (2) and equations (Al) through (A8) are utilized in 
the main program. Input-output is, of course, accomplished here. 

Subroutine ITER (MS2501) is used to compute the points of intersection of 
the ray with the window surfaces. Equations (4) and (5) and (A9) through (A14) 


6 














are used in this subroutine. The calling statement is CALL ITER (X, Y, SM, 
Cl, K, DELZ, DELTAP) where: 

X x coordinate 

Y y coordinate 

SM SM(I, J, K) = Mi matrices for math models 

Cl Cl (I) 1 = 1, 3 = I x , Iy, I z components of I 

K index 

DELZ AZ 

DELTAP scale factor 

SM and Cl are arrays. X, Y, SM, Cl, K, and DELTAP are input from the main 
program. DELZ is computed in ITER and returned. New values of X and Y 
are computed also and returned. 

Subroutine NORMAL (MS2502) is used to compute the normals to the window 
surfaces at the points of intersection. Equations (6), ( 7 ), and (8) are used 
in this subroutine. The calling statement is CALL NORMAL (X, Y, SM, CN, K, 
DELTAP) where: 

X x coordinate 

Y y coordinate 

SM SM(I, J, K) = Mi matrices for the math models 

CN CN(I) = N x , Ny, N z components of N 

K index 

DELTAP scale factor 

SM and CN are arrays. The components CN (I ) of N are computed in the 
subroutine and returned to main program. Others are input from the main 
program. 

Subroutine REFRC (MS2303) uses equation (1) to compute the refracted ray 
in each medium through which the ray passes. The calling statement is CALL 
REFRC (Cl, CN, QRI , CR) where 


Cl 

Cl (I) = Ix, Iy, I z 

components 

of 

I 

CN 

CN(I) = N x , N y , N z 

components 

of 

N 
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n k+l 

QRI ratio of refractive indices 

% 

CR CR(I) = R x , Ry, R z components of R 

Cl, CN, and CR are arrays. The components CR(I) of the refracted ray R 

are computed in this routine and the other parameters come from the main 

program. 

Subroutine XPYXM (MS2504) is a matrix multiplying subroutine used in 
subroutines ITER and NORMAL. The calling sequence is (A, B, C, NRA, NCA, 

NRB, NCB, J) where A and B are the matrices to be multiplied and NRA, NCA, 

NRB, and NCB are, respectively, the number of rows in A, columns in A, rows 

in B, and columns in B; C is the result of the matrix multiplication. The 
index J determines the multiplication to be performed. For 

J = 1 C = AB 

J = 2 C = AB t 

j = 3 C = A T B 

Listings and flow charts for the main program and subroutines appear in 

appendixes B and C. 

; PROGRAM USAGE 


The program is written in the FORTRAN IV computer language. It operates 
at Ames Research Center on an IBM 7094 computer under the IBJOB Processor of 
the IBSYS operating system, version 13. The program requires 2752 8 or 1514 10 
storage locations exclusive of the FORTRAN system subroutines required. 


DECK MAKEUP 


$JOB 

CARD 


$ IBJOB 

CARD 


MS2300 

DECK 

Main program 

MS2301 

DECK 

ITER 

MS2302 

DECK 

NORMAL 

MS2303 ' 

DECK 

REFRC 

MS2304 

DECK 

XPYXM 

$DATA 

CARD 


INPUT 

CARDS 



PROGRAM INPUTS 

Inputs to the ray trace program can be divided into three categories: 
constants, math model matrices, and individual ray data. The constants 
describe the window configuration and its environment. The math model 
matrices describe the window surface shapes. The individual ray data describe 
the rays to be traced through the window. The data are all input in card 
form. The input cards are described below. 
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Constants 


Columns 


Symbol 

Definition 

Card 1 

Format 

110 


1-10 


N 

Number of window surfaces 

Card 2 

Format 

8E10.0 


1-10 


D(l) 

Distances from reference planes 

11-20 


D (2) 

to other reference planes 

21-30 


D (3) 


31-40 


D(4) 


41-50 


D(5) 


Card 3 

Format 

8E10.0 


1-10 


DELTAP (1) 

Scale factors at interfaces 

11-20 


DELTAP (2) 

between air and glass 

21-30 


DELTAP (3) 


31-40 


DELTAP (4) 


41-50 


DELTAP (5) 


51-60 


DELTAP (6) 


Card 4 

Format 

8E10.0 


1-10 


Ri(l) 

Refractive indices of each 

11-20 


RI(2) 

medium through which ray is 

21-30 


Ri (3) 

traced 

31-40 


RI(4) 


41-50 


Ri (5) 


51-60 


RI (6) 


61-70 


RI (7) 



See figure 3. 


RP[ - RP2 RP3 '■ RP4 RP5 RP6 



Figure 3.- Definition of reference planes, 
refractive indices, and input distances. 
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Math Model Matrices 


Cards 5-9 Columns of math model matrix 1 

Format 5E15.8 J = 1, 5 

Columns 1-15 Mj 4 

16-30 M 2 a 

31-45 M 3 i 

46-60 M 4 ’^ 

61-75 M’i 
5 > J 

Cards 10-14 Columns of math model matrix 2 
if N = 2 Same format 

Cards 15-19 Columns of math model matrix 3 
if N = 3 Same format 

Cards 20-24 Columns of math model matrix 4 
if N = 4 Same format 

Cards 25-29 Columns of math model matrix 5 
if N = 5 Same format 

Cards 30-34 Columns of math model matrix 6 
if N = 6 Same format 


Individual Ray Data 

Cards 35-°° Format 8E10.0 

1-10 ALPHA I | Ray orientation angles a^, 6-j. 

11-20 DELTA I J 

21-30 E(l)) x, y, z coordinates of 

31-40 E(2) i reference position from which 

41-50 E(3) ) ray originates 

There will be one card for each separate ray to be traced for a given set 
of constants and surface matrices. When one ray is traced through the window, 
the program looks for another to trace and the program is terminated when 
there are no more rays to trace. As presently written, the program considers 
only one set of constants and surface matrices. If it is desired to consider 
another set of constants or surface matrices or both, the program must be 
rerun. The input cards as described above are for a three-pane window. If 
the window has more or less than three panes, there would be a corresponding 
increase or decrease in the number of cards. 
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PROGRAM OUTPUTS 


All the program inputs are printed out in the same sequence that they are 
input. In addition, for each input ray the following quantities are output: 


Symbol 


ALPHAR 

DELTAR 

DELALP 

DELDEL 

DELINC 

x 

y 

z 


a r | orientation angles of final refracted ray on 
(5 r j outside of window 

a i ~ a r = Aa 
6^ - = A6 

Angle between incident ray and final refracted ray 

Xf, yp, Zf coordinates of point where refracted ray 
leaves outermost window surface 


A sample data printout is included with this writeup. 


The program requires 0.2 to 0.3 second of computation time to complete a 
single trace through a three-pane window. 


TAPES 


Logical tape 5 is used for input and logical tape 6 for output. 

SAMPLE CASE DESCRIPTION 


The printout for a sample ray trace follows this section. In the sample 
case the window consisted of one pane of glass 0.3 inch thick with a pressure 
of 5 psi imposed on it from within the spacecraft. The window is circular 
with a 7-inch diameter. Two rays, both with = 0° and 6^ = 60°, were 
traced through the window. One ray intersected the inner window surface at 
point X = 0, Y = 0 and the other at point X = 6, Y = 0. 
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CONCLUDING REMARKS 


A FORTRAN IV computer program has been written for tracing rays 
computationally through spacecraft windows. 

Program equations are written in vector-matrix form for ease in computa- 
tion. The program requires 0.2 to 0.3 second to trace a ray through a three- 
pane window. The program determines the angular deviations of the ray as it 
passes through the window. 

Rays may be traced through any window whose surface shapes can be 
described mathematically in polynomial form. 

The ray trace program has been used extensively to trace rays through 
Gemini spacecraft optical windows as discussed in reference 1. The program is 
currently being used to trace rays through Apollo windows and also through 
generalized windows of various sizes and shapes. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, Dec. 5, 1969 
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APPENDIX A 


COMPUTATION OF THE POINT OF INTERSECTION OF A RAY WITH 
THE REFERENCE PLANE AND THE WINDOW SURFACE 


This appendix describes the determination of the intersection of the 
assumed interior ray I with the first reference plane RP1 and also the 


determination of the intersection of 



I with the first window surface. This 
is an iterative procedure, and is illus- 
trated in figure 4. Here the window 
coordinate system is given as (ijk) ; 
the vector from the origin locating 
either the eye of the observer, or some 
reference point on the instrument is 
given as E, some initial ray vector as 


Figure 4.- Geometry of the scheme for iterating 
from the reference plane to the window 
surface. 


and from figure 4 


where 


Then from the third equation 


Ti , which terminates in the 

x-y plane, 

and a vector in the x-y plane closing 
with Ii as d. Let 

A A A 

E - E x i + Eyj + E z k 

(Al) 

ii = ai + gj + yk 

(A2) 

Ti = ill 

(A3) 

where Ii is unknown. Now 
the form 

A -A 

d is of 

+ y x j + Ok 

(A4) 

E + III 

E x + Il° \ 

(A5) 

Ey + i x e > 
E z - I lY ) 

(A6) 

II 

1 

(A7) 


and from the first two 
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X 1 - E x 


y l = E y 

Now, if we set 

AIj. = 0 (A9) 

compute 

z i = ApY^tMjX! (A10) 

S\ 

and project the increment of z\ above the Ij vector onto the I direction 

A = (zi - A I i I • k)k • I = (zi - AIiy)y (All) 

Now form 

AI 2 = Ali + A (A12) 

and compute 

x 2 = xi + AI • i = xi + Act (A13) 

y 2 = y i + A * ' j = y l + ( A14 ) 

These values of x 2 and y 2 are then used in the equation for z (eq. (A10)), 
and the computation is continued until j A { < e. At present e = l.OxlO" 6 in 
the program. At this time the last values computed for x, and y, and z 
are the coordinates of the point of intersection of our ray with the refract- 
ing surface, to the accuracy we have specified for e. The basic scheme of 
this iteration is that by successively projecting the value of (z - AI • k) 
onto the vector I, we approach the point of intersection with the window 
surface. This system appeals to be stable for all continuous smooth surface 
functions . 
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APPENDIX B 


PROGRAM LISTING 
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